{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# **ICESat-2 Active Subglacial Lakes in Antarctica**\n",
    "\n",
    "Finding subglacial lakes that are draining or filling under the ice!\n",
    "They can be detected with ICESat-2 data, as significant changes in height\n",
    "(> 1 metre) over a relatively short duration (< 1 year), i.e. a high rate of\n",
    "elevation change over time (dhdt).\n",
    "\n",
    "In this notebook, we'll use some neat tools to help us examine the lakes:\n",
    "- To find active subglacial lake boundaries,\n",
    "use an *unsupervised clustering* technique\n",
    "- To see ice surface elevation trends at a higher temporal resolution (< 3 months),\n",
    "perform *crossover track error analysis* on intersecting ICESat-2 tracks\n",
    "\n",
    "To speed up analysis on millions of points,\n",
    "we will use state of the art GPU algorithms enabled by RAPIDS AI libraries,\n",
    "or parallelize the processing across our HPC's many CPU cores using Dask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "import subprocess\n",
    "\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\"\n",
    "\n",
    "import cudf\n",
    "import cuml\n",
    "import dask\n",
    "import dask.array\n",
    "import geopandas as gpd\n",
    "import hvplot.xarray\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import panel as pn\n",
    "import pygmt\n",
    "import scipy.spatial\n",
    "import shapely.geometry\n",
    "import tqdm\n",
    "import xarray as xr\n",
    "import zarr\n",
    "\n",
    "import deepicedrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tag: str = \"X2SYS\"\n",
    "os.environ[\"X2SYS_HOME\"] = os.path.abspath(tag)\n",
    "client = dask.distributed.Client(\n",
    "    n_workers=64, threads_per_worker=1, env={\"X2SYS_HOME\": os.environ[\"X2SYS_HOME\"]}\n",
    ")\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "min_date, max_date = (\"2018-10-14\", \"2020-07-16\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(\"ATLXI/df_dhdt_antarctica.parquet\"):\n",
    "    zarrarray = zarr.open_consolidated(store=f\"ATLXI/ds_dhdt_antarctica.zarr\", mode=\"r\")\n",
    "    _ = deepicedrain.ndarray_to_parquet(\n",
    "        ndarray=zarrarray,\n",
    "        parquetpath=\"ATLXI/df_dhdt_antarctica.parquet\",\n",
    "        variables=[\"x\", \"y\", \"dhdt_slope\", \"referencegroundtrack\", \"h_corr\"],\n",
    "        dropnacols=[\"dhdt_slope\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Read in Antarctic Drainage Basin Boundaries shapefile into a GeoDataFrame\n",
    "ice_boundaries: gpd.GeoDataFrame = (\n",
    "    deepicedrain.catalog.measures_antarctic_boundaries.read()\n",
    ")\n",
    "drainage_basins: gpd.GeoDataFrame = ice_boundaries.query(expr=\"TYPE == 'GR'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in ICESat-2 data (x, y, dhdt) and do initial trimming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trimmed 212404341 -> 53002087\n"
     ]
    }
   ],
   "source": [
    "# Read in raw x, y, dhdt_slope and referencegroundtrack data into the GPU\n",
    "cudf_raw: cudf.DataFrame = cudf.read_parquet(\n",
    "    filepath_or_buffer=\"ATLXI/df_dhdt_antarctica.parquet\",\n",
    "    columns=[\"x\", \"y\", \"dhdt_slope\", \"referencegroundtrack\"],\n",
    ")\n",
    "# Filter to points with dhdt that is less than -0.2 m/yr or more than +0.2 m/yr\n",
    "cudf_many = cudf_raw.loc[abs(cudf_raw.dhdt_slope) > 0.2]\n",
    "print(f\"Trimmed {len(cudf_raw)} -> {len(cudf_many)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clip outlier values to 3 sigma (standard deviations) from mean\n",
    "_mean = cudf_many.dhdt_slope.mean()\n",
    "_std = cudf_many.dhdt_slope.std()\n",
    "cudf_many.dhdt_slope.clip(\n",
    "    lower=np.float32(_mean - 3 * _std), upper=np.float32(_mean + 3 * _std), inplace=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Label ICESat-2 points according to their drainage basin\n",
    "\n",
    "Uses Point in Polygon.\n",
    "For each point, find out which Antarctic Drainage Basin they are in.\n",
    "This will also remove the points on floating (FR) ice shelves and islands (IS),\n",
    "so that we keep only points on the grounded (GR) ice regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trimmed 53002087 -> 27848833\n"
     ]
    }
   ],
   "source": [
    "# Use point in polygon to label points according to the drainage basins they fall in\n",
    "cudf_many[\"drainage_basin\"]: cudf.Series = deepicedrain.point_in_polygon_gpu(\n",
    "    points_df=cudf_many, poly_df=drainage_basins\n",
    ")\n",
    "X_many = cudf_many.dropna()  # drop points that are not in a drainage basin\n",
    "print(f\"Trimmed {len(cudf_many)} -> {len(X_many)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# Find Active Subglacial Lake clusters\n",
    "\n",
    "Uses Density-based spatial clustering of applications with noise (DBSCAN)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subglacial Lake Finder algorithm\n",
    "\n",
    "For each Antarctic drainage basin:\n",
    "\n",
    "1. Select all points with significant elevation change over time (dhdt)\n",
    "  - Specifically, the (absolute) dhdt value should be\n",
    "    2x the median (absolute) dhdt for that drainage basin\n",
    "  - E.g. if median dhdt for basin is 0.35 m/yr,\n",
    "    we choose points that have dhdt > 0.70 m/yr\n",
    "2. Run unsupervised clustering to pick out active subglacial lakes\n",
    "  - Split into draining (-dhdt) and filling (+dhdt) points first\n",
    "  - Use DBSCAN algorithm to cluster points into groups,\n",
    "    with an eps (distance) of 3 km and minimum sample size of 250 points\n",
    "3. Check each potential point cluster to see if it meets active lake criteria\n",
    "  1. Build a convex hull 'lake' polygon around clustered points\n",
    "  2. Check that the 'lake' has significant elevation change relative to outside\n",
    "    - For the area in the 5 km buffer region **outside** the 'lake' polygon:\n",
    "       - Find median dhdt (outer_dhdt)\n",
    "       - Find median absolute deviation of dhdt values (outer_mad)\n",
    "    - For the area **inside** the 'lake' polygon:\n",
    "       - Find median dhdt (inner_dhdt)\n",
    "    - If the potential lake shows an elevation change that is more than\n",
    "      3x the surrounding deviation of background elevation change,\n",
    "      we infer that this is likely an active subglacial 'lake'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/198 [00:00<01:34,  2.07it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "75298 rows at Academy above ± 0.55 m/yr\n",
      "1 draining and 6 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 21%|██        | 41/198 [00:04<00:34,  4.56it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "79410 rows at Totten above ± 0.77 m/yr\n",
      "3 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 28%|██▊       | 56/198 [00:05<00:13, 10.85it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33313 rows at Cook above ± 0.58 m/yr\n",
      "0 draining and 1 filling lakes found\n",
      "22517 rows at David above ± 0.58 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 29%|██▉       | 58/198 [00:07<00:48,  2.90it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "107148 rows at Mercer above ± 0.66 m/yr\n",
      "5 draining and 10 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 32%|███▏      | 63/198 [00:17<03:42,  1.65s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "317073 rows at Pine_Island above ± 0.86 m/yr\n",
      "2 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 34%|███▍      | 67/198 [00:21<03:30,  1.60s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "243245 rows at Thwaites above ± 0.72 m/yr\n",
      "3 draining and 6 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 36%|███▌      | 71/198 [00:24<02:29,  1.17s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "144101 rows at Whillans above ± 0.66 m/yr\n",
      "5 draining and 16 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 36%|███▋      | 72/198 [00:24<02:04,  1.01it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36739 rows at Kamb above ± 0.63 m/yr\n",
      "2 draining and 6 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 37%|███▋      | 73/198 [00:24<01:34,  1.32it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10402 rows at Leverett above ± 0.81 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 37%|███▋      | 74/198 [00:26<02:10,  1.05s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96475 rows at Scott above ± 0.65 m/yr\n",
      "3 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 39%|███▉      | 77/198 [00:27<01:24,  1.43it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51012 rows at Amundsen above ± 0.69 m/yr\n",
      "3 draining and 0 filling lakes found\n",
      "25640 rows at Shackleton above ± 0.76 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 40%|████      | 80/198 [00:28<00:50,  2.32it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76667 rows at Beardmore above ± 0.69 m/yr\n",
      "2 draining and 1 filling lakes found\n",
      "22900 rows at Lennox-King above ± 0.77 m/yr\n",
      "2 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 41%|████▏     | 82/198 [00:28<00:41,  2.78it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50217 rows at Nimrod above ± 0.64 m/yr\n",
      "3 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 43%|████▎     | 85/198 [00:29<00:33,  3.36it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67558 rows at Byrd above ± 0.56 m/yr\n",
      "7 draining and 3 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 45%|████▌     | 90/198 [00:30<00:27,  3.86it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48095 rows at MacAyeal above ± 0.60 m/yr\n",
      "5 draining and 2 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 47%|████▋     | 93/198 [00:30<00:21,  4.99it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51268 rows at Rutford above ± 0.76 m/yr\n",
      "0 draining and 2 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 48%|████▊     | 96/198 [00:31<00:20,  5.08it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "86860 rows at Evans above ± 0.70 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 55%|█████▌    | 109/198 [00:31<00:08, 11.11it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45012 rows at Slessor above ± 0.59 m/yr\n",
      "11 draining and 4 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 57%|█████▋    | 112/198 [00:33<00:18,  4.78it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34026 rows at Support_Force above ± 0.56 m/yr\n",
      "3 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 58%|█████▊    | 114/198 [00:34<00:24,  3.44it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "94847 rows at Foundation above ± 0.55 m/yr\n",
      "1 draining and 5 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 60%|██████    | 119/198 [00:34<00:14,  5.28it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16586 rows at Lambert above ± 0.52 m/yr\n",
      "0 draining and 1 filling lakes found\n",
      "16511 rows at Mellor above ± 0.53 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 67%|██████▋   | 133/198 [00:35<00:05, 12.52it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "38850 rows at Institute above ± 0.56 m/yr\n",
      "5 draining and 2 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 70%|██████▉   | 138/198 [00:36<00:11,  5.30it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30250 rows at Kohler above ± 1.79 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 71%|███████   | 140/198 [00:37<00:10,  5.77it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30867 rows at Ross_East2 above ± 0.96 m/yr\n",
      "3 draining and 0 filling lakes found\n",
      "37890 rows at Bowman_Strom_Live_Axel-Heigerg above ± 0.99 m/yr\n",
      "4 draining and 1 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 72%|███████▏  | 142/198 [00:37<00:13,  4.20it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45785 rows at Sulzberger above ± 0.89 m/yr\n",
      "2 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 73%|███████▎  | 144/198 [00:38<00:10,  5.00it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5264 rows at Lucchitta_Velasco above ± 0.94 m/yr\n",
      "1 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 74%|███████▎  | 146/198 [00:40<00:26,  1.99it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "144411 rows at Getz above ± 1.29 m/yr\n",
      "2 draining and 0 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 83%|████████▎ | 164/198 [00:42<00:06,  5.27it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "61080 rows at Recovery above ± 0.53 m/yr\n",
      "4 draining and 4 filling lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 198/198 [00:44<00:00,  4.48it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total of 159 subglacial lakes found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Subglacial lake finder\n",
    "activelakes: dict = {\n",
    "    \"basin_name\": [],  # Antarctic drainage basin name\n",
    "    \"num_points\": [],  # Number of clustered data points\n",
    "    \"outer_dhdt\": [],  # Median elevation change over time (dhdt) outside of lake\n",
    "    \"outer_std\": [],  # Standard deviation of dhdt outside of lake\n",
    "    \"outer_mad\": [],  # Median absolute deviation of dhdt outside of lake\n",
    "    \"inner_dhdt\": [],  # Median elev change over time (dhdt) inside of lake bounds\n",
    "    \"maxabsdhdt\": [],  # Maximum absolute dhdt value inside of lake boundary\n",
    "    \"refgtracks\": [],  # Pipe-delimited list of ICESat-2 reference ground tracks\n",
    "    \"geometry\": [],  # Shapely Polygon geometry holding lake boundary coordinates\n",
    "}\n",
    "basin_name: str = \"Whillans\"  # Set a basin name here\n",
    "basins = drainage_basins[drainage_basins.NAME == basin_name].index  # one specific basin\n",
    "basins: pd.core.indexes.numeric.Int64Index = drainage_basins.index  # run on all basins\n",
    "for basin_index in tqdm.tqdm(iterable=basins):\n",
    "    # Initial data cleaning, filter to rows that are in the drainage basin\n",
    "    basin = drainage_basins.loc[basin_index]\n",
    "    X_local = X_many.loc[X_many.drainage_basin == basin.NAME]  # .reset_index(drop=True)\n",
    "\n",
    "    # Get points with dhdt_slope higher than 2x the median dhdt_slope for the basin\n",
    "    # E.g. if median dhdt_slope is 0.35 m/yr, then we cluster points over 0.70 m/yr\n",
    "    abs_dhdt = X_local.dhdt_slope.abs()\n",
    "    tolerance: float = 2 * abs_dhdt.median()\n",
    "    X = X_local.loc[abs_dhdt > tolerance]\n",
    "\n",
    "    if len(X) <= 1000:  # don't run on too few points\n",
    "        continue\n",
    "\n",
    "    # Run unsupervised clustering separately on draining and filling lakes\n",
    "    # Draining lake points have negative labels (e.g. -1, -2, 3),\n",
    "    # Filling lake points have positive labels (e.g. 1, 2, 3),\n",
    "    # Noise points have NaN labels (i.e. NaN)\n",
    "    cluster_vars = [\"x\", \"y\", \"dhdt_slope\"]\n",
    "    draining_lake_labels = -deepicedrain.find_clusters(\n",
    "        X=X.loc[X.dhdt_slope < 0][cluster_vars]\n",
    "    )\n",
    "    filling_lake_labels = deepicedrain.find_clusters(\n",
    "        X=X.loc[X.dhdt_slope > 0][cluster_vars]\n",
    "    )\n",
    "    lake_labels = cudf.concat(objs=[draining_lake_labels, filling_lake_labels])\n",
    "    lake_labels: cudf.Series = lake_labels.sort_index()\n",
    "\n",
    "    # Checking all potential subglacial lakes in a basin\n",
    "    clusters: cudf.Series = lake_labels.unique()\n",
    "    for cluster_label in clusters.to_array():\n",
    "        # Store attribute and geometry information of each active lake\n",
    "        lake_points: cudf.DataFrame = X.loc[lake_labels == cluster_label]\n",
    "\n",
    "        # More data cleaning, dropping clusters with too few points\n",
    "        try:\n",
    "            assert len(lake_points) > 100\n",
    "        except AssertionError:\n",
    "            lake_labels = lake_labels.replace(to_replace=cluster_label, value=None)\n",
    "            continue\n",
    "\n",
    "        multipoint: shapely.geometry.MultiPoint = shapely.geometry.MultiPoint(\n",
    "            points=lake_points[[\"x\", \"y\"]].as_matrix()\n",
    "        )\n",
    "        convexhull: shapely.geometry.Polygon = multipoint.convex_hull\n",
    "\n",
    "        # Filter out (most) false positive subglacial lakes\n",
    "        # Check that elevation change over time in lake is anomalous to outside\n",
    "        # The 5000 m distance from lake boundary setting is empirically based on\n",
    "        # Smith et al. 2009's methodology at https://doi.org/10.3189/002214309789470879\n",
    "        outer_ring_buffer = convexhull.buffer(distance=5000) - convexhull\n",
    "        X_local[\"in_donut_ring\"] = deepicedrain.point_in_polygon_gpu(\n",
    "            points_df=X_local,\n",
    "            poly_df=gpd.GeoDataFrame({\"name\": True, \"geometry\": [outer_ring_buffer]}),\n",
    "        )\n",
    "        outer_points = X_local.dropna(subset=\"in_donut_ring\")\n",
    "        outer_dhdt: float = outer_points.dhdt_slope.median()\n",
    "\n",
    "        outer_std: float = outer_points.dhdt_slope.std()\n",
    "        outer_mad: float = scipy.stats.median_abs_deviation(\n",
    "            x=outer_points.dhdt_slope.to_pandas()\n",
    "        )\n",
    "\n",
    "        inner_dhdt: float = lake_points.dhdt_slope.median()\n",
    "        X_local = X_local.drop(labels=\"in_donut_ring\", axis=\"columns\")\n",
    "\n",
    "        # If lake interior's median dhdt value is within 3 median absolute deviations\n",
    "        # of the lake exterior's dhdt value, we remove the lake label\n",
    "        # I.e. skip if above background change not significant enough\n",
    "        # Inspired by Kim et al. 2016's methodology at https://doi.org/10.5194/tc-10-2971-2016\n",
    "        if abs(inner_dhdt - outer_dhdt) < 3 * outer_mad:\n",
    "            lake_labels = lake_labels.replace(to_replace=cluster_label, value=None)\n",
    "            continue\n",
    "\n",
    "        maxabsdhdt: float = (\n",
    "            lake_points.dhdt_slope.max()\n",
    "            if cluster_label > 0  # positive label = filling\n",
    "            else lake_points.dhdt_slope.min()  # negative label = draining\n",
    "        )\n",
    "        refgtracks: str = \"|\".join(\n",
    "            map(str, lake_points.referencegroundtrack.unique().to_pandas())\n",
    "        )\n",
    "\n",
    "        # Save key variables to dictionary that will later go into geodataframe\n",
    "        activelakes[\"basin_name\"].append(basin.NAME)\n",
    "        activelakes[\"num_points\"].append(len(lake_points))\n",
    "        activelakes[\"outer_dhdt\"].append(outer_dhdt)\n",
    "        activelakes[\"outer_std\"].append(outer_std)\n",
    "        activelakes[\"outer_mad\"].append(outer_mad)\n",
    "        activelakes[\"inner_dhdt\"].append(inner_dhdt)\n",
    "        activelakes[\"maxabsdhdt\"].append(maxabsdhdt)\n",
    "        activelakes[\"refgtracks\"].append(refgtracks)\n",
    "        activelakes[\"geometry\"].append(convexhull)\n",
    "\n",
    "    # Calculate total number of lakes found for one drainage basin\n",
    "    clusters: cudf.Series = lake_labels.unique()\n",
    "    n_draining, n_filling = (clusters < 0).sum(), (clusters > 0).sum()\n",
    "    if n_draining + n_filling > 0:\n",
    "        print(f\"{len(X)} rows at {basin.NAME} above ± {tolerance:.2f} m/yr\")\n",
    "        print(f\"{n_draining} draining and {n_filling} filling lakes found\")\n",
    "\n",
    "if len(activelakes[\"geometry\"]) >= 1:\n",
    "    gdf = gpd.GeoDataFrame(activelakes, crs=\"EPSG:3031\")\n",
    "    basename = \"antarctic_subglacial_lakes\"  # f\"temp_{basin_name.lower()}_lakes\"  #\n",
    "    gdf.to_file(filename=f\"{basename}_3031.geojson\", driver=\"GeoJSON\")\n",
    "    gdf.to_crs(crs={\"init\": \"epsg:4326\"}).to_file(\n",
    "        filename=f\"{basename}_4326.geojson\", driver=\"GeoJSON\"\n",
    "    )\n",
    "\n",
    "print(f\"Total of {len(gdf)} subglacial lakes found\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize lakes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Concatenate XY points with labels, and move data from GPU to CPU\n",
    "X: cudf.DataFrame = cudf.concat(objs=[X, lake_labels], axis=\"columns\")\n",
    "X_ = X.to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Plot clusters on a map in colour, noise points/outliers as small dots\n",
    "fig = pygmt.Figure()\n",
    "n_clusters_ = len(X_.cluster_id.unique()) - 1  # No. of clusters minus noise (NaN)\n",
    "sizes = (X_.cluster_id.isna()).map(arg={True: 0.01, False: 0.1})\n",
    "pygmt.makecpt(cmap=\"polar\", series=(-1, 3, 2), color_model=\"+cDrain,Fill\", reverse=True)\n",
    "fig.plot(\n",
    "    x=X_.x,\n",
    "    y=X_.y,\n",
    "    sizes=sizes,\n",
    "    style=\"cc\",\n",
    "    color=pd.cut(x=X_.cluster_id, bins=(-np.inf, 0, np.inf), labels=[-1, 1]),\n",
    "    cmap=True,\n",
    "    frame=[\n",
    "        f'WSne+t\"Estimated number of lake clusters at {basin.NAME}: {n_clusters_}\"',\n",
    "        'xafg+l\"Polar Stereographic X (m)\"',\n",
    "        'yafg+l\"Polar Stereographic Y (m)\"',\n",
    "    ],\n",
    ")\n",
    "basinx, basiny = basin.geometry.exterior.coords.xy\n",
    "fig.plot(x=basinx, y=basiny, pen=\"thinnest,-\")\n",
    "fig.colorbar(position='JMR+w2c/0.5c+e+m+n\"Unclassified\"', L=\"i0.5c\")\n",
    "fig.savefig(fname=f\"figures/subglacial_lake_clusters_at_{basin.NAME}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Select a subglacial lake to examine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Save or load dhdt data from Parquet file\n",
    "placename: str = \"siple_coast\"  # \"Recovery\"  # \"Whillans\"\n",
    "try:\n",
    "    drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys=\"NAME\")\n",
    "    region: deepicedrain.Region = deepicedrain.Region.from_gdf(\n",
    "        gdf=drainage_basins.loc[placename], name=\"Recovery Basin\"\n",
    "    )\n",
    "except KeyError:\n",
    "    pass\n",
    "df_dhdt: cudf.DataFrame = cudf.read_parquet(\n",
    "    f\"ATLXI/df_dhdt_{placename.lower()}.parquet\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Antarctic subglacial lake polygons with EPSG:3031 coordinates\n",
    "antarctic_lakes: gpd.GeoDataFrame = deepicedrain.catalog.subglacial_lakes.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "basin_name                                             Whillans\n",
      "num_points                                                 3276\n",
      "outer_dhdt                                             0.377611\n",
      "outer_std                                               1.42712\n",
      "outer_mad                                              0.083278\n",
      "inner_dhdt                                              2.12607\n",
      "maxabsdhdt                                              8.65662\n",
      "refgtracks    74|135|196|266|327|388|516|577|638|769|830|101...\n",
      "geometry      POLYGON ((-444681.0351994165 -545210.454471163...\n",
      "Name: 50, dtype: object\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"300\" height=\"300\" viewBox=\"-456991.1908163784 -545964.7144366237 20365.019067435118 18110.406576156733\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,-1073819.0222970909)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"135.76679378290078\" opacity=\"0.6\" d=\"M -444681.0351994165,-545210.4544711632 L -445059.4273992341,-545139.2030781403 L -447613.01038991945,-544185.9610993187 L -448604.6076458619,-543785.8544671827 L -448890.77904354926,-543669.5833568621 L -451986.0162835539,-542018.9649696343 L -454210.39452298917,-540418.8503920811 L -454408.2814011552,-540246.0395154255 L -455385.85207539203,-539162.6159726225 L -456224.89981433837,-537787.3535748171 L -456236.9308509178,-537500.1188677208 L -454509.0218480694,-532882.2389798661 L -454459.2975782437,-532851.9066176108 L -448142.3565408698,-529956.985603917 L -441807.6895563269,-528608.5678259275 L -441356.8655731373,-528717.5355566369 L -439391.13700329117,-530087.2909355733 L -438887.0564588413,-530681.5730842056 L -438879.5685567435,-530699.8250218468 L -437919.08106515603,-534731.9868097329 L -437380.43171440385,-542757.777147059 L -437400.2001665644,-542812.5670977917 L -437460.8738999268,-542976.4327690548 L -437481.0986066531,-543031.0538226315 L -438289.55547334306,-544026.4356033215 L -438338.1695958391,-544058.5196375257 L -444681.0351994165,-545210.4544711632 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f2190a02c10>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Choose one draining/filling lake\n",
    "draining: bool = False\n",
    "placename: str = \"Whillans\"  # \"Slessor\"  # \"Kamb\"  # \"Mercer\"  #\n",
    "lakes: gpd.GeoDataFrame = antarctic_lakes.query(expr=\"basin_name == @placename\")\n",
    "lake = lakes.loc[lakes.inner_dhdt.idxmin() if draining else lakes.inner_dhdt.idxmax()]\n",
    "lake = lakes.query(expr=\"inner_dhdt < 0\" if draining else \"inner_dhdt > 0\").loc[48]\n",
    "lakedict = {\n",
    "    21: \"Mercer 2b\",  # filling lake\n",
    "    40: \"Lower Subglacial Lake Conway\",  # draining lake\n",
    "    41: \"Subglacial Lake Conway\",  # draining lake\n",
    "    48: \"Subglacial Lake Whillans\",  # filling lake\n",
    "    50: \"Whillans IX\",  # filling lake\n",
    "    63: \"Kamb 1\",  # filling lake\n",
    "    65: \"Kamb 12\",  # filling lake\n",
    "    97: \"MacAyeal 1\",  # draining lake\n",
    "    109: \"Slessor 45\",  # draining lake\n",
    "    116: \"Slessor 23\",  # filling lake\n",
    "    151: \"Recovery IV\",  # draining lake\n",
    "    156: \"Recovery 2\",  # filling lake\n",
    "}\n",
    "region = deepicedrain.Region.from_gdf(gdf=lake, name=lakedict[lake.name])\n",
    "\n",
    "print(lake)\n",
    "lake.geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Subset data to lake of interest\n",
    "placename: str = region.name.lower().replace(\" \", \"_\")\n",
    "df_lake: cudf.DataFrame = region.subset(data=df_dhdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an interpolated ice surface elevation grid for each ICESat-2 cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 6/6 [00:03<00:00,  1.71it/s]\n"
     ]
    }
   ],
   "source": [
    "# Generate gridded time-series of ice elevation over lake\n",
    "cycles: tuple = (3, 4, 5, 6, 7, 8)\n",
    "os.makedirs(name=f\"figures/{placename}\", exist_ok=True)\n",
    "ds_lake: xr.Dataset = deepicedrain.spatiotemporal_cube(\n",
    "    table=df_lake.to_pandas(),\n",
    "    placename=placename,\n",
    "    cycles=cycles,\n",
    "    folder=f\"figures/{placename}\",\n",
    ")\n",
    "ds_lake.to_netcdf(path=f\"figures/{placename}/xyht_{placename}.nc\", mode=\"w\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elevation limits are: (273.2418518066406, 310.9927062988281)\n"
     ]
    }
   ],
   "source": [
    "# Get 3D grid_region (xmin/xmax/ymin/ymax/zmin/zmax),\n",
    "# and calculate normalized z-values as Elevation delta relative to Cycle 3\n",
    "z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max()))  # original z limits\n",
    "grid_region: tuple = region.bounds() + z_limits\n",
    "\n",
    "ds_lake_diff: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z\n",
    "z_diff_limits: tuple = (float(ds_lake_diff.z.min()), float(ds_lake_diff.z.max()))\n",
    "diff_grid_region: np.ndarray = np.append(arr=grid_region[:4], values=z_diff_limits)\n",
    "\n",
    "print(f\"Elevation limits are: {z_limits}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D plots of gridded ice surface elevation over time\n",
    "for cycle in tqdm.tqdm(iterable=cycles):\n",
    "    time_nsec: pd.Timestamp = df_lake[f\"utc_time_{cycle}\"].to_pandas().mean()\n",
    "    time_sec: str = np.datetime_as_string(arr=time_nsec.to_datetime64(), unit=\"s\")\n",
    "\n",
    "    grid = f\"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc\"\n",
    "    points = pd.DataFrame(\n",
    "        data=np.vstack(lake.geometry.boundary.coords.xy).T, columns=(\"x\", \"y\")\n",
    "    )\n",
    "    outline_points = pygmt.grdtrack(points=points, grid=grid, newcolname=\"z\")\n",
    "    outline_points[\"z\"] = outline_points.z.fillna(value=outline_points.z.median())\n",
    "\n",
    "    fig = deepicedrain.plot_icesurface(\n",
    "        grid=grid,  # ds_lake.sel(cycle_number=cycle).z\n",
    "        grid_region=grid_region,\n",
    "        diff_grid=ds_lake_diff.sel(cycle_number=cycle).z,\n",
    "        diff_grid_region=diff_grid_region,\n",
    "        track_points=df_lake[[\"x\", \"y\", f\"h_corr_{cycle}\"]].dropna().as_matrix(),\n",
    "        outline_points=outline_points,\n",
    "        azimuth=157.5,  # 202.5  # 270\n",
    "        elevation=45,  # 60\n",
    "        title=f\"{region.name} at Cycle {cycle} ({time_sec})\",\n",
    "    )\n",
    "    fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Make an animated GIF of changing ice surface from the PNG files\n",
    "# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname}\n",
    "gif_fname: str = (\n",
    "    f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif\"\n",
    ")\n",
    "subprocess.check_call(\n",
    "    [\n",
    "        \"convert\",\n",
    "        \"-delay\",\n",
    "        \"120\",\n",
    "        \"-loop\",\n",
    "        \"0\",\n",
    "        f\"figures/{placename}/dsm_*.png\",\n",
    "        gif_fname,\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HvPlot 2D interactive view of ice surface elevation grids over each ICESat-2 cycle\n",
    "dashboard: pn.layout.Column = pn.Column(\n",
    "    ds_lake.hvplot.image(x=\"x\", y=\"y\", clim=z_limits, cmap=\"gist_earth\", data_aspect=1)\n",
    "    # * ds_lake.hvplot.contour(x=\"x\", y=\"y\", clim=z_limits, data_aspect=1)\n",
    ")\n",
    "dashboard.show(port=30227)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Along track plots of ice surface elevation change over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a few Reference Ground tracks to look at\n",
    "rgts: list = [int(rgt) for rgt in lake.refgtracks.split(\"|\")]\n",
    "print(f\"Looking at Reference Ground Tracks: {rgts}\")\n",
    "os.makedirs(name=f\"figures/{placename}\", exist_ok=True)\n",
    "\n",
    "track_dict: dict = {}\n",
    "rgt_groups = df_lake.groupby(by=\"referencegroundtrack\")\n",
    "for rgt, df_rgt_wide in tqdm.tqdm(rgt_groups, total=len(rgt_groups.groups.keys())):\n",
    "    df_rgt: pd.DataFrame = deepicedrain.wide_to_long(\n",
    "        df=df_rgt_wide.to_pandas(), stubnames=[\"h_corr\", \"utc_time\"], j=\"cycle_number\"\n",
    "    )\n",
    "\n",
    "    # Split one referencegroundtrack into 3 laser pair tracks pt1, pt2, pt3\n",
    "    df_rgt[\"pairtrack\"]: pd.Series = pd.cut(\n",
    "        x=df_rgt.y_atc, bins=[-np.inf, -100, 100, np.inf], labels=(\"pt1\", \"pt2\", \"pt3\")\n",
    "    )\n",
    "    pt_groups = df_rgt.groupby(by=\"pairtrack\")\n",
    "    for pairtrack, df_ in pt_groups:\n",
    "        if len(df_) > 0:\n",
    "            rgtpair = f\"{rgt:04d}_{pairtrack}\"\n",
    "            track_dict[rgtpair] = df_\n",
    "\n",
    "            # Transect plot along a reference ground track\n",
    "            fig = deepicedrain.plot_alongtrack(\n",
    "                df=df_, rgtpair=rgtpair, regionname=region.name, oldtonew=draining\n",
    "            )\n",
    "            fig.savefig(\n",
    "                fname=f\"figures/{placename}/alongtrack_{placename}_{rgtpair}.png\"\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# Crossover Track Analysis\n",
    "\n",
    "To increase the temporal resolution of\n",
    "our ice elevation change analysis\n",
    "(i.e. at time periods less than\n",
    "the 91 day repeat cycle of ICESat-2),\n",
    "we can look at the locations where the\n",
    "ICESat-2 tracks intersect and get the\n",
    "height values there!\n",
    "Uses [pygmt.x2sys_cross](https://www.pygmt.org/v0.2.0/api/generated/pygmt.x2sys_cross.html).\n",
    "\n",
    "References:\n",
    "- Wessel, P. (2010). Tools for analyzing intersecting tracks: The x2sys package.\n",
    "Computers & Geosciences, 36(3), 348–354. https://doi.org/10.1016/j.cageo.2009.05.009"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize X2SYS database in the X2SYS/ICESAT2 folder\n",
    "pygmt.x2sys_init(\n",
    "    tag=\"ICESAT2\",\n",
    "    fmtfile=f\"{tag}/ICESAT2/xyht\",\n",
    "    suffix=\"tsv\",\n",
    "    units=[\"de\", \"se\"],  # distance in metres, speed in metres per second\n",
    "    gap=\"d250e\",  # distance gap up to 250 metres allowed\n",
    "    force=True,\n",
    "    verbose=\"q\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Run crossover analysis on all tracks\n",
    "rgts, tracks = track_dict.keys(), track_dict.values()\n",
    "# Parallelized paired crossover analysis\n",
    "futures: list = []\n",
    "for rgt1, rgt2 in itertools.combinations(rgts, r=2):\n",
    "    # skip if same referencegroundtrack but different laser pair\n",
    "    # as they are parallel and won't cross\n",
    "    if rgt1[:4] == rgt2[:4]:\n",
    "        continue\n",
    "    track1 = track_dict[rgt1][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
    "    track2 = track_dict[rgt2][[\"x\", \"y\", \"h_corr\", \"utc_time\"]]\n",
    "    future = client.submit(\n",
    "        key=f\"{rgt1}x{rgt2}\",\n",
    "        func=pygmt.x2sys_cross,\n",
    "        tracks=[track1, track2],\n",
    "        tag=\"ICESAT2\",\n",
    "        # region=[-460000, -400000, -560000, -500000],\n",
    "        interpolation=\"l\",  # linear interpolation\n",
    "        coe=\"e\",  # external crossovers\n",
    "        trackvalues=True,  # Get track 1 height (h_1) and track 2 height (h_2)\n",
    "        # trackvalues=False,  # Get crossover error (h_X) and mean height value (h_M)\n",
    "        # outfile=\"xover_236_562.tsv\"\n",
    "    )\n",
    "    futures.append(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "crossovers: dict = {}\n",
    "for f in tqdm.tqdm(\n",
    "    iterable=dask.distributed.as_completed(futures=futures), total=len(futures)\n",
    "):\n",
    "    if f.status != \"error\":  # skip those track pairs which don't intersect\n",
    "        crossovers[f.key] = f.result().dropna().reset_index(drop=True)\n",
    "\n",
    "df_cross: pd.DataFrame = pd.concat(objs=crossovers, names=[\"track1_track2\", \"id\"])\n",
    "df: pd.DataFrame = df_cross.reset_index(level=\"track1_track2\").reset_index(drop=True)\n",
    "# Report on how many unique crossover intersections there were\n",
    "# df.plot.scatter(x=\"x\", y=\"y\")  # quick plot of our crossover points\n",
    "print(\n",
    "    f\"{len(df.groupby(by=['x', 'y']))} crossover intersection point locations found \"\n",
    "    f\"with {len(df)} crossover height-time pairs \"\n",
    "    f\"over {len(tracks)} tracks\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate crossover error\n",
    "df[\"h_X\"]: pd.Series = df.h_2 - df.h_1  # crossover error (i.e. height difference)\n",
    "df[\"t_D\"]: pd.Series = df.t_2 - df.t_1  # elapsed time in ns (i.e. time difference)\n",
    "ns_in_yr: int = 365.25 * 24 * 60 * 60 * 1_000_000_000  # nanoseconds in a year\n",
    "df[\"dhdt\"]: pd.Series = df.h_X / (df.t_D.astype(np.int64) / ns_in_yr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Get some summary statistics of our crossover errors\n",
    "sumstats: pd.DataFrame = df[[\"h_X\", \"t_D\", \"dhdt\"]].describe()\n",
    "# Find location with highest absolute crossover error, and most sudden height change\n",
    "max_h_X: pd.Series = df.iloc[np.nanargmax(df.h_X.abs())]  # highest crossover error\n",
    "max_dhdt: pd.Series = df.iloc[df.dhdt.argmax()]  # most sudden change in height"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D Map view of crossover points\n",
    "\n",
    "Bird's eye view of the crossover points\n",
    "overlaid on top of the ICESat-2 tracks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# 2D plot of crossover locations\n",
    "var: str = \"h_X\"\n",
    "fig = pygmt.Figure()\n",
    "# Setup basemap\n",
    "plotregion = pygmt.info(table=df[[\"x\", \"y\"]], spacing=1000)\n",
    "pygmt.makecpt(cmap=\"batlow\", series=[sumstats[var][\"25%\"], sumstats[var][\"75%\"]])\n",
    "# Map frame in metre units\n",
    "fig.basemap(frame=\"f\", region=plotregion, projection=\"X8c\")\n",
    "# Plot actual track points in green\n",
    "for track in tracks:\n",
    "    tracklabel = f\"{track.iloc[0].referencegroundtrack} {track.iloc[0].pairtrack}\"\n",
    "    fig.plot(\n",
    "        x=track.x,\n",
    "        y=track.y,\n",
    "        pen=\"thinnest,green,.\",\n",
    "        style=f'qN+1:+l\"{tracklabel}\"+f3p,Helvetica,darkgreen',\n",
    "    )\n",
    "# Plot crossover point locations\n",
    "fig.plot(x=df.x, y=df.y, color=df.h_X, cmap=True, style=\"c0.1c\", pen=\"thinnest\")\n",
    "# Plot lake boundary in blue\n",
    "lakex, lakey = lake.geometry.exterior.coords.xy\n",
    "fig.plot(x=lakex, y=lakey, pen=\"thin,blue,-.\")\n",
    "# Map frame in kilometre units\n",
    "fig.basemap(\n",
    "    frame=[\n",
    "        f'WSne+t\"Crossover points at {region.name}\"',\n",
    "        'xaf+l\"Polar Stereographic X (km)\"',\n",
    "        'yaf+l\"Polar Stereographic Y (km)\"',\n",
    "    ],\n",
    "    region=plotregion / 1000,\n",
    "    projection=\"X8c\",\n",
    ")\n",
    "fig.colorbar(position=\"JMR\", frame=['x+l\"Crossover Error\"', \"y+lm\"])\n",
    "fig.savefig(f\"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Crossover Elevation time-series\n",
    "\n",
    "Plot elevation change over time at:\n",
    "\n",
    "1. One single crossover point location\n",
    "2. Many crossover locations over an area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tidy up dataframe first using pd.wide_to_long\n",
    "# I.e. convert 't_1', 't_2', 'h_1', 'h_2' columns into just 't' and 'h'.\n",
    "df_th: pd.DataFrame = deepicedrain.wide_to_long(\n",
    "    df=df.loc[:, [\"track1_track2\", \"x\", \"y\", \"t_1\", \"t_2\", \"h_1\", \"h_2\"]],\n",
    "    stubnames=[\"t\", \"h\"],\n",
    "    j=\"track\",\n",
    ")\n",
    "df_th = df_th.drop_duplicates(ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot at single location with **maximum** absolute crossover height error (max_h_X)\n",
    "df_max = df_th.query(expr=\"x == @max_h_X.x & y == @max_h_X.y\").sort_values(by=\"t\")\n",
    "track1, track2 = df_max.track1_track2.iloc[0].split(\"x\")\n",
    "print(f\"{round(max_h_X.h_X, 2)} metres height change at {max_h_X.x}, {max_h_X.y}\")\n",
    "plotregion = np.array(\n",
    "    [df_max.t.min(), df_max.t.max(), *pygmt.info(table=df_max[[\"h\"]], spacing=2.5)[:2]]\n",
    ")\n",
    "plotregion += np.array([-pd.Timedelta(2, unit=\"W\"), +pd.Timedelta(2, unit=\"W\"), 0, 0])\n",
    "\n",
    "fig = pygmt.Figure()\n",
    "with pygmt.config(\n",
    "    FONT_ANNOT_PRIMARY=\"9p\", FORMAT_TIME_PRIMARY_MAP=\"abbreviated\", FORMAT_DATE_MAP=\"o\"\n",
    "):\n",
    "    fig.basemap(\n",
    "        projection=\"X12c/8c\",\n",
    "        region=plotregion,\n",
    "        frame=[\n",
    "            f'WSne+t\"Max elevation change over time at {region.name}\"',\n",
    "            \"pxa1Of1o+lDate\",  # primary time axis, 1 mOnth annotation and minor axis\n",
    "            \"sx1Y\",  # secondary time axis, 1 Year intervals\n",
    "            'yaf+l\"Elevation at crossover (m)\"',\n",
    "        ],\n",
    "    )\n",
    "fig.text(\n",
    "    text=f\"Track {track1} and {track2} crossover\",\n",
    "    position=\"TC\",\n",
    "    offset=\"jTC0c/0.2c\",\n",
    "    V=\"q\",\n",
    ")\n",
    "# Plot data points\n",
    "fig.plot(x=df_max.t, y=df_max.h, style=\"c0.15c\", color=\"darkblue\", pen=\"thin\")\n",
    "# Plot dashed line connecting points\n",
    "fig.plot(x=df_max.t, y=df_max.h, pen=f\"faint,blue,-\")\n",
    "fig.savefig(\n",
    "    f\"figures/{placename}/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png\"\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all crossover height points over time over the lake area\n",
    "fig = deepicedrain.plot_crossovers(df=df_th, regionname=region.name)\n",
    "fig.savefig(f\"figures/{placename}/crossover_many_{placename}_{min_date}_{max_date}.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all crossover height points over time over the lake area\n",
    "# with height values normalized to 0 from the first observation date\n",
    "normfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()\n",
    "df_th[\"h_norm\"] = df_th.groupby(by=\"track1_track2\").h.transform(func=normfunc)\n",
    "\n",
    "fig = deepicedrain.plot_crossovers(\n",
    "    df=df_th,\n",
    "    regionname=region.name,\n",
    "    elev_var=\"h_norm\",\n",
    "    elev_filter=3 * abs(df.h_X).median(),\n",
    ")\n",
    "fig.savefig(\n",
    "    f\"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png\"\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "formats": "ipynb,py:hydrogen"
  },
  "kernelspec": {
   "display_name": "deepicedrain",
   "language": "python",
   "name": "deepicedrain"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
